home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / minpack / hybrj.f < prev    next >
Text File  |  1996-07-19  |  14KB  |  442 lines
  1.       SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,MODE,
  2.      *                 FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,WA1,WA2,
  3.      *                 WA3,WA4)
  4.       INTEGER N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR
  5.       DOUBLE PRECISION XTOL,FACTOR
  6.       DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N),DIAG(N),R(LR),
  7.      *                 QTF(N),WA1(N),WA2(N),WA3(N),WA4(N)
  8.       EXTERNAL FCN
  9. C     **********
  10. C
  11. C     SUBROUTINE HYBRJ
  12. C
  13. C     THE PURPOSE OF HYBRJ IS TO FIND A ZERO OF A SYSTEM OF
  14. C     N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
  15. C     OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A
  16. C     SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN.
  17. C
  18. C     THE SUBROUTINE STATEMENT IS
  19. C
  20. C       SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,
  21. C                        MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,
  22. C                        WA1,WA2,WA3,WA4)
  23. C
  24. C     WHERE
  25. C
  26. C       FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
  27. C         CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST
  28. C         BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER
  29. C         CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
  30. C
  31. C         SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
  32. C         INTEGER N,LDFJAC,IFLAG
  33. C         DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N)
  34. C         ----------
  35. C         IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND
  36. C         RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC.
  37. C         IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND
  38. C         RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC.
  39. C         ---------
  40. C         RETURN
  41. C         END
  42. C
  43. C         THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
  44. C         THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ.
  45. C         IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
  46. C
  47. C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  48. C         OF FUNCTIONS AND VARIABLES.
  49. C
  50. C       X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
  51. C         AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
  52. C         CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
  53. C
  54. C       FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
  55. C         THE FUNCTIONS EVALUATED AT THE OUTPUT X.
  56. C
  57. C       FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
  58. C         ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION
  59. C         OF THE FINAL APPROXIMATE JACOBIAN.
  60. C
  61. C       LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
  62. C         WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
  63. C
  64. C       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION
  65. C         OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE
  66. C         ITERATES IS AT MOST XTOL.
  67. C
  68. C       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
  69. C         OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1
  70. C         HAS REACHED MAXFEV.
  71. C
  72. C       DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE
  73. C         BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG
  74. C         MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS IMPLICIT
  75. C         (MULTIPLICATIVE) SCALE FACTORS FOR THE VARIABLES.
  76. C
  77. C       MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE
  78. C         VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,
  79. C         THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER
  80. C         VALUES OF MODE ARE EQUIVALENT TO MODE = 1.
  81. C
  82. C       FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE
  83. C         INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF
  84. C         FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE
  85. C         TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE
  86. C         INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.
  87. C
  88. C       NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED
  89. C         PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,
  90. C         FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST
  91. C         ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND
  92. C         IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE
  93. C         FOR PRINTING. FVEC AND FJAC SHOULD NOT BE ALTERED.
  94. C         IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS OF FCN
  95. C         WITH IFLAG = 0 ARE MADE.
  96. C
  97. C       INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
  98. C         TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
  99. C         VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
  100. C         INFO IS SET AS FOLLOWS.
  101. C
  102. C         INFO = 0   IMPROPER INPUT PARAMETERS.
  103. C
  104. C         INFO = 1   RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES
  105. C                    IS AT MOST TOL.
  106. C
  107. C         INFO = 2   NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS
  108. C                    REACHED MAXFEV.
  109. C
  110. C         INFO = 3   XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
  111. C                    THE APPROXIMATE SOLUTION X IS POSSIBLE.
  112. C
  113. C         INFO = 4   ITERATION IS NOT MAKING GOOD PROGRESS, AS
  114. C                    MEASURED BY THE IMPROVEMENT FROM THE LAST
  115. C                    FIVE JACOBIAN EVALUATIONS.
  116. C
  117. C         INFO = 5   ITERATION IS NOT MAKING GOOD PROGRESS, AS
  118. C                    MEASURED BY THE IMPROVEMENT FROM THE LAST
  119. C                    TEN ITERATIONS.
  120. C
  121. C       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
  122. C         CALLS TO FCN WITH IFLAG = 1.
  123. C
  124. C       NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
  125. C         CALLS TO FCN WITH IFLAG = 2.
  126. C
  127. C       R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE
  128. C         UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION
  129. C         OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE.
  130. C
  131. C       LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
  132. C         (N*(N+1))/2.
  133. C
  134. C       QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
  135. C         THE VECTOR (Q TRANSPOSE)*FVEC.
  136. C
  137. C       WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N.
  138. C
  139. C     SUBPROGRAMS CALLED
  140. C
  141. C       USER-SUPPLIED ...... FCN
  142. C
  143. C       MINPACK-SUPPLIED ... DOGLEG,DPMPAR,ENORM,
  144. C                            QFORM,QRFAC,R1MPYQ,R1UPDT
  145. C
  146. C       FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,MOD
  147. C
  148. C     MINPACK. VERSION OF SEPTEMBER 1979.
  149. C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
  150. C
  151. C     **********
  152.       INTEGER I,IFLAG,ITER,J,JM1,L,NCFAIL,NCSUC,NSLOW1,NSLOW2
  153.       INTEGER IWA(1)
  154.       LOGICAL JEVAL,SING
  155.       DOUBLE PRECISION ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,
  156.      *                 PRERED,P1,P5,P001,P0001,RATIO,SUM,TEMP,XNORM,
  157.      *                 ZERO
  158.       DOUBLE PRECISION DPMPAR,ENORM
  159.       DATA ONE,P1,P5,P001,P0001,ZERO
  160.      *     /1.0D0,1.0D-1,5.0D-1,1.0D-3,1.0D-4,0.0D0/
  161. C
  162. C     EPSMCH IS THE MACHINE PRECISION.
  163. C
  164.       EPSMCH = DPMPAR(1)
  165. C
  166.       INFO = 0
  167.       IFLAG = 0
  168.       NFEV = 0
  169.       NJEV = 0
  170. C
  171. C     CHECK THE INPUT PARAMETERS FOR ERRORS.
  172. C
  173.       IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. XTOL .LT. ZERO
  174.      *    .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO
  175.      *    .OR. LR .LT. (N*(N + 1))/2) GO TO 300
  176.       IF (MODE .NE. 2) GO TO 20
  177.       DO 10 J = 1, N
  178.          IF (DIAG(J) .LE. ZERO) GO TO 300
  179.    10    CONTINUE
  180.    20 CONTINUE
  181. C
  182. C     EVALUATE THE FUNCTION AT THE STARTING POINT
  183. C     AND CALCULATE ITS NORM.
  184. C
  185.       IFLAG = 1
  186.       CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
  187.       NFEV = 1
  188.       IF (IFLAG .LT. 0) GO TO 300
  189.       FNORM = ENORM(N,FVEC)
  190. C
  191. C     INITIALIZE ITERATION COUNTER AND MONITORS.
  192. C
  193.       ITER = 1
  194.       NCSUC = 0
  195.       NCFAIL = 0
  196.       NSLOW1 = 0
  197.       NSLOW2 = 0
  198. C
  199. C     BEGINNING OF THE OUTER LOOP.
  200. C
  201.    30 CONTINUE
  202.          JEVAL = .TRUE.
  203. C
  204. C        CALCULATE THE JACOBIAN MATRIX.
  205. C
  206.          IFLAG = 2
  207.          CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
  208.          NJEV = NJEV + 1
  209.          IF (IFLAG .LT. 0) GO TO 300
  210. C
  211. C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
  212. C
  213.          CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)
  214. C
  215. C        ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
  216. C        TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
  217. C
  218.          IF (ITER .NE. 1) GO TO 70
  219.          IF (MODE .EQ. 2) GO TO 50
  220.          DO 40 J = 1, N
  221.             DIAG(J) = WA2(J)
  222.             IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
  223.    40       CONTINUE
  224.    50    CONTINUE
  225. C
  226. C        ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
  227. C        AND INITIALIZE THE STEP BOUND DELTA.
  228. C
  229.          DO 60 J = 1, N
  230.             WA3(J) = DIAG(J)*X(J)
  231.    60       CONTINUE
  232.          XNORM = ENORM(N,WA3)
  233.          DELTA = FACTOR*XNORM
  234.          IF (DELTA .EQ. ZERO) DELTA = FACTOR
  235.    70    CONTINUE
  236. C
  237. C        FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.
  238. C
  239.          DO 80 I = 1, N
  240.             QTF(I) = FVEC(I)
  241.    80       CONTINUE
  242.          DO 120 J = 1, N
  243.             IF (FJAC(J,J) .EQ. ZERO) GO TO 110
  244.             SUM = ZERO
  245.             DO 90 I = J, N
  246.                SUM = SUM + FJAC(I,J)*QTF(I)
  247.    90          CONTINUE
  248.             TEMP = -SUM/FJAC(J,J)
  249.             DO 100 I = J, N
  250.                QTF(I) = QTF(I) + FJAC(I,J)*TEMP
  251.   100          CONTINUE
  252.   110       CONTINUE
  253.   120       CONTINUE
  254. C
  255. C        COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.
  256. C
  257.          SING = .FALSE.
  258.          DO 150 J = 1, N
  259.             L = J
  260.             JM1 = J - 1
  261.             IF (JM1 .LT. 1) GO TO 140
  262.             DO 130 I = 1, JM1
  263.                R(L) = FJAC(I,J)
  264.                L = L + N - I
  265.   130          CONTINUE
  266.   140       CONTINUE
  267.             R(L) = WA1(J)
  268.             IF (WA1(J) .EQ. ZERO) SING = .TRUE.
  269.   150       CONTINUE
  270. C
  271. C        ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.
  272. C
  273.          CALL QFORM(N,N,FJAC,LDFJAC,WA1)
  274. C
  275. C        RESCALE IF NECESSARY.
  276. C
  277.          IF (MODE .EQ. 2) GO TO 170
  278.          DO 160 J = 1, N
  279.             DIAG(J) = DMAX1(DIAG(J),WA2(J))
  280.   160       CONTINUE
  281.   170    CONTINUE
  282. C
  283. C        BEGINNING OF THE INNER LOOP.
  284. C
  285.   180    CONTINUE
  286. C
  287. C           IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
  288. C
  289.             IF (NPRINT .LE. 0) GO TO 190
  290.             IFLAG = 0
  291.             IF (MOD(ITER-1,NPRINT) .EQ. 0)
  292.      *         CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
  293.             IF (IFLAG .LT. 0) GO TO 300
  294.   190       CONTINUE
  295. C
  296. C           DETERMINE THE DIRECTION P.
  297. C
  298.             CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)
  299. C
  300. C           STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
  301. C
  302.             DO 200 J = 1, N
  303.                WA1(J) = -WA1(J)
  304.                WA2(J) = X(J) + WA1(J)
  305.                WA3(J) = DIAG(J)*WA1(J)
  306.   200          CONTINUE
  307.             PNORM = ENORM(N,WA3)
  308. C
  309. C           ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
  310. C
  311.             IF (ITER .EQ. 1) DELTA = DMIN1(DELTA,PNORM)
  312. C
  313. C           EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
  314. C
  315.             IFLAG = 1
  316.             CALL FCN(N,WA2,WA4,FJAC,LDFJAC,IFLAG)
  317.             NFEV = NFEV + 1
  318.             IF (IFLAG .LT. 0) GO TO 300
  319.             FNORM1 = ENORM(N,WA4)
  320. C
  321. C           COMPUTE THE SCALED ACTUAL REDUCTION.
  322. C
  323.             ACTRED = -ONE
  324.             IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
  325. C
  326. C           COMPUTE THE SCALED PREDICTED REDUCTION.
  327. C
  328.             L = 1
  329.             DO 220 I = 1, N
  330.                SUM = ZERO
  331.                DO 210 J = I, N
  332.                   SUM = SUM + R(L)*WA1(J)
  333.                   L = L + 1
  334.   210             CONTINUE
  335.                WA3(I) = QTF(I) + SUM
  336.   220          CONTINUE
  337.             TEMP = ENORM(N,WA3)
  338.             PRERED = ZERO
  339.             IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2
  340. C
  341. C           COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
  342. C           REDUCTION.
  343. C
  344.             RATIO = ZERO
  345.             IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED
  346. C
  347. C           UPDATE THE STEP BOUND.
  348. C
  349.             IF (RATIO .GE. P1) GO TO 230
  350.                NCSUC = 0
  351.                NCFAIL = NCFAIL + 1
  352.                DELTA = P5*DELTA
  353.                GO TO 240
  354.   230       CONTINUE
  355.                NCFAIL = 0
  356.                NCSUC = NCSUC + 1
  357.                IF (RATIO .GE. P5 .OR. NCSUC .GT. 1)
  358.      *            DELTA = DMAX1(DELTA,PNORM/P5)
  359.                IF (DABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5
  360.   240       CONTINUE
  361. C
  362. C           TEST FOR SUCCESSFUL ITERATION.
  363. C
  364.             IF (RATIO .LT. P0001) GO TO 260
  365. C
  366. C           SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
  367. C
  368.             DO 250 J = 1, N
  369.                X(J) = WA2(J)
  370.                WA2(J) = DIAG(J)*X(J)
  371.                FVEC(J) = WA4(J)
  372.   250          CONTINUE
  373.             XNORM = ENORM(N,WA2)
  374.             FNORM = FNORM1
  375.             ITER = ITER + 1
  376.   260       CONTINUE
  377. C
  378. C           DETERMINE THE PROGRESS OF THE ITERATION.
  379. C
  380.             NSLOW1 = NSLOW1 + 1
  381.             IF (ACTRED .GE. P001) NSLOW1 = 0
  382.             IF (JEVAL) NSLOW2 = NSLOW2 + 1
  383.             IF (ACTRED .GE. P1) NSLOW2 = 0
  384. C
  385. C           TEST FOR CONVERGENCE.
  386. C
  387.             IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1
  388.             IF (INFO .NE. 0) GO TO 300
  389. C
  390. C           TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
  391. C
  392.             IF (NFEV .GE. MAXFEV) INFO = 2
  393.             IF (P1*DMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3
  394.             IF (NSLOW2 .EQ. 5) INFO = 4
  395.             IF (NSLOW1 .EQ. 10) INFO = 5
  396.             IF (INFO .NE. 0) GO TO 300
  397. C
  398. C           CRITERION FOR RECALCULATING JACOBIAN.
  399. C
  400.             IF (NCFAIL .EQ. 2) GO TO 290
  401. C
  402. C           CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN
  403. C           AND UPDATE QTF IF NECESSARY.
  404. C
  405.             DO 280 J = 1, N
  406.                SUM = ZERO
  407.                DO 270 I = 1, N
  408.                   SUM = SUM + FJAC(I,J)*WA4(I)
  409.   270             CONTINUE
  410.                WA2(J) = (SUM - WA3(J))/PNORM
  411.                WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM)
  412.                IF (RATIO .GE. P0001) QTF(J) = SUM
  413.   280          CONTINUE
  414. C
  415. C           COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.
  416. C
  417.             CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)
  418.             CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)
  419.             CALL R1MPYQ(1,N,QTF,1,WA2,WA3)
  420. C
  421. C           END OF THE INNER LOOP.
  422. C
  423.             JEVAL = .FALSE.
  424.             GO TO 180
  425.   290    CONTINUE
  426. C
  427. C        END OF THE OUTER LOOP.
  428. C
  429.          GO TO 30
  430.   300 CONTINUE
  431. C
  432. C     TERMINATION, EITHER NORMAL OR USER IMPOSED.
  433. C
  434.       IF (IFLAG .LT. 0) INFO = IFLAG
  435.       IFLAG = 0
  436.       IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
  437.       RETURN
  438. C
  439. C     LAST CARD OF SUBROUTINE HYBRJ.
  440. C
  441.       END
  442.